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ABSTRACT 

The linear normal-mode stratorotational instability (SRI) is analytically reexamined in the in- 
viscid limit where the length scales of horizontal disturbances are large compared their vertical 
and radial counterparts. Boundary conditions different than channel walls are also considered. 
This quasi-hydrostatic, semi-geostrophic (QHSG) approximation allows one to examine the 
effect of a vertically varying Brunt- Vaisaila frequency, N 2 . It is found that the normal-mode 
instability persists when iV 2 increases quadratically with respect to the disc vertical coordi- 
nate. However we also find that the SRI seems to exist in this inviscid QHSG extreme only 
for channel wall conditions: when one or both of the reflecting walls are removed there is no 
instability in the asymptotic limit explored here. It is also found that only exponential-type 
SRI modes (as defined by Dubrulle et al. 2005) exist under these conditions. These equations 
also admit non-normal mode behaviour. Fixed Lagrangian pressure conditions on both radial 
boundaries predicts there to be no normal mode behaviour in the QHSG limit. The mathemat- 
ical relationship between the results obtained here and that of the classic Eady (1949) problem 
for baroclinic instability is drawn. We conjecture as to the mathematical/physical nature of the 
SRI. 

The general linear problem, analyzed without approximation in the context of the Boussi- 
nesq equations, admits a potential vorticity-like quantity that is advectively conserved by the 
shear. Its existence means that a continuous spectrum is a generic feature of this system. It also 
implies that in places where the Brunt- Vaisaila frequency becomes dominant the linearized 
flow may two-dimensionalize by advectively conserving its vertical vorticity. 
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The question whether or not hydrodynamic activity can 
emerge in protoplanetary discs is experiencing a Renaissance. 
Given the absence of an inflection point in the basic Keplerian flow 
of discs, it was natural for investigations, beginning in the early 
90's, to consider other physical effects as a possible source of su- 
percritical linear instabilities. The MRI (Balbus, 2003) an instabil- 
ity (non-conservative) involving the joint interplay of rotation and 
magnetic effects, has proven itself to be a viable linear mechanism 
which could lead to globally sustained activity in discs. 

However, fresh analysis of simplified models of protoplan- 
etary discs, like the shearing sheet approximation (Goldreich & 
Lynden-Bell, 1965) utilized in many hydrodynamic and magneto- 
hydrodynamic investigations of circumstellar Keplerian discs, have 
shown that a number of alternative routes to long-term activity, 
both linear and nonlinear, can occur for purely hydrodynamic dis- 
turbances. Bracco et al (1999) demonstrated that anticy clonic vor- 
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tices can live for long periods of time before ultimately decaying 
away in viscous two-dimensional, global disc models of incom- 
pressible flow executing Keplerian motion. This work has since led 
to many other investigations which consider an array of dynami- 
cal processes, sometimes even linearly stable, that could operate in 
a purely hydrodynamic manner and these include (but not limited 
to): large amplitude defects of vorticity in the flow (Li et al. 2000) 
leading to inflection point "secondary" - instabilities, steady forc- 
ing of 2D linearized global viscous disc disturbances which have 
shown to lead to strong transiently growing structures and patterns 
(Iounnou & Kakouris, 2002). Transient growth in 2D disturbances, 
an effect primarily driven by the Orr-mechanism of vortex-tilting 
(Orr 1907) and which is fundamentally a non-normal mode effect 
(Schmid & Henningson, 2001, Chagelishvilli et al. 2003, Tevzadze 
et al., 2003, Yecko 2004) can lead to sustained nonlinear activity in 
2D, including robust and coherent anticyclones (Umurhan & Regev 
2004). 

Recently, Barranco and Marcus (2005), considering an anal- 
ogous set of equations describing the shearing- sheet limit of a cir- 
cumstellar disc, showed that nonlinear disturbances in 3D, where 
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the disc vertical direction and effect of gravity are included, lead 
to long-lived vortex structures with vorticity pointing parallel to 
the direction of gravity: something akin to steady anticyclones per- 
sisting away from the disc midplane. This is a striking feature be- 
cause it implies that in these regimes the flow behaves nearly two- 
dimensionally. This goes against implications of the shearing sheet 
simulations of Balbus, Hawley & Stone (1996) and Hawley, Balbus 
& Winters (1999) which have suggested that there is slim possibil- 
ity of sustained dynamics in purely hydrodynamic systems subject 
to strong rotation and shear. 

The effect of buoyancy on disturbances has only recently been 
reconsidered. According to the Solberg-Hoiland criterion it seemed 
that a sufficient condition for a naturally occurring non-magnetic 
instability in a disc is for there to be an adverse entropy gradient 
in the direction normal to the disc plane. \ However the work of 
Yavneh et al. (2002) has recently turned this wisdom on its head. 
They demonstrated that a linear stability analysis of the classic Tay- 
lor Couette problem (that is, flow between two concentric rotating 
cylinders) in which the fluid is stratified (in the direction parallel 
to the cylinder axes) leads to normal-mode instability even if the 
stratification is stable to buoyant oscillations. 

Dubrulle et al., 2005 (D05) showed that 3D linearized normal- 
mode disturbances of stratified local sections of a protoplanetary 
disc can manifest this instability. The effect, termed the Stratoro- 
tational Instability (SRI), has been proposed as a mechanism, free 
of non-conservative effects like magnetic field interactions, which 
could actively operate in the disc, possibly lead to turbulence 
(D05), and has been suggested to be relevant in connection to the 
steady vortices observed in the simulations of Barranco & Marcus 
(2005) (Shalybkov & Rudiger, 2005). 

The SRI effect, as revealed by Yavneh et al (2000), Shalybkov 
& Rudiger (2005) and in D05, was considered for flows in which 
the radial boundaries of the system are rigid: in other words, un- 
der boundary conditions in which there is no normal flow across 
the bounding channel walls. § In order to affect a tractable analysis 
the investigations thus far have assumed that the disc vertical com- 
ponent of gravity, gz, and the vertical gradient of the basic state 
temperature, d z Tb (or, equivalently, the basic state entropy gradi- 
ent, d z Sb) are constant with respect to the disc vertical coordinate. 
Although this may be a good approximation for atmospheric flows 
like the Earth's (Pedlosky, 1987), this is not an accurate represen- 
tation of conditions near the midplane of a circumstellar disc. 

Furthermore, the question as to what are reasonable radial 
boundary conditions for shearing-sheet sections of Keplerian discs 
remains very open. Shearing-sheet sections of Keplerian discs, cen- 
tered about the disc midplane, have both a g and a d z Tt (d z Sb), 
whose product is proportional to the square of the Brunt-Vaisaila 
frequency, N 2 , which varies quadratically with respect to the height 
from the midplane, z. It is natural then to ask (i) whether or not the 
SRI is a dynamical effect driven primarily by the reflecting property 
of the system's radial boundaries and, (ii) whether or not it persists 
under conditions in which the Brunt-Vaisaila frequency varies with 
respect to the vertical coordinate - something one would expect in 
a protoplanetary disc. 

A purpose of this study is to revisit, primarily via asymptotic 
means, the properties of the SRI instability as discussed by D05. 

X The acoustic instability of Papaloizou & Pringle (1984,1985) is not con- 
sidered "natural" in this sense because it relies on the presence of artificial 
boundary conditions. See Section 5 for further discussion on this matter. 
§ Note, however, that D05 also consider stress-free and periodic boundaries 
in the radial directions. 



The main motivation of this inquiry is to see whether or not the 
no-flow radial boundary conditions employed by their study and/or 
whether or not the assumption of the constancy of N 2 significantly 
alters the SRI effect uncovered by previous investigations. 

We would like to make it clear here that the inclusion of strat- 
ification in quantities like the vertical component of gravity and 
the temperature/entropy gradients mathematically results in a nor- 
mal mode problem who is primarily non-separable in the radial and 
vertical coordinates. It is for this reason a general analysis of this 
physical scenario is quite difficult for both analytical and numerical 
reasons. 

We approach these questions by reconsidering the SRI within 
the context of two model equations describing a local section of a 
circumstellar disc. We will experiment with the effects of differing 
boundary conditions. 

In Section 1, we introduce the first of the equations, the Large- 
Shearing Box (LSB), which are the basic equations appropriate to 
shearing-sheet sections of accretion discs centered about their mid- 
planes. These are then used to motivate a second, simpler, model 
equation set: the incompressible Boussinesq equations (BE) in in- 
viscid three dimensional rotating plane Couette flow (or rpCf, see 
Yecko, 2004). This simpler set is mathematically equivalent to the 
set studied, in the inviscid limit, in D05. Boundary conditions are 
considered which share a common property in the BE, namely 
those that cause the total disturbance energy, called E, to change 
solely due to the energy exchanged between disturbances and the 
shear through the Reynolds stress term and not due to work done 
on the system from outside. This is achieved by enforcing periodic- 
ity in the azimuthal direction: either periodicity in the disc vertical 
direction (when constant N 2 is assumed) or zero normal (vertical) 
velocity fluctuations at some fiducial vertical disc boundary (when 
N varies with disc height). There a number of radial boundary 
conditions which achieves these objective and we investigate these: 
(a) no-normal flow at the radial boundaries, sometimes referred to 
in this manuscript as "channel-wall conditions", (b) no Lagrangian 
variation of the pressure on the moving radial boundaries, (c) or 
some mixed combination of these two. 

In Section 2 normal mode solutions of the BE are considered. 
First it is shown that there exists an advectively conserved quan- 
tity of the flow. The conservation of this quantity in certain limits 
implies that the linearized flow is nearly two-dimensional - in that 
the vertical vorticity is the dominant quantity that is preserved by 
the shear adverted flow. We then proceed towards a normal mode 
analysis by initially assuming the constancy of the Brunt-Vaisaila 
frequency, N 2 . We asymptotically analyze disturbances whose ra- 
dial (x) and vertical length scales are dwarfed by comparison to the 
azimuthal scales (y) - a circumstance which is refereed to as the 
quasi-hydrostatic semi-geostrophic approximation (QHSG). This 
limit predisposes the resulting equations into admitting simple an- 
alytical solutions. This asymptotic limit shares many of the quali- 
ties of the quasi- WKB analysis considered in D05. Implementing 
the variety of boundary conditions enumerated above shows that, 
in this limit, only the no normal (radial) flow boundary conditions 
admit unstable normal modes. We also find that there is always 
a continuous spectrum in this problem, irrespective of the radial 
boundary conditions, and we briefly discuss some of its features. 
The section is rounded out by relaxing the constancy of N 2 and 
considering the case where it varies quadratically with respect to 
the disc height in this same QHSG limit. There appears to be a 
potential vorticity-like quantity which is advectively conserved in 
the linear limit. Moreover, it turns out that it is possible to con- 
struct separable solutions (in x and z) for the normal mode prob- 
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lem. The main conclusion from this is that the normal mode sta- 
bility behaviour is unaffected by the stratification of TV 2 . Unlike 
in D05, where the SRI appears for both "exponentiaF'-modes and 
"oscillatory"-modes, the results of this section show that the SRI 
only occurs for "exponentiaF'-modes. 

In Section 3 the same QHSG approximation is applied directly 
to the LSB. These yield, similarly, a conserved potential vortic- 
ity like quantity but which now includes the effects of weak com- 
pressibility and a finite soundspeed. ^[ The equations for the normal- 
modes are also separable here and it is found, again, that the sta- 
bility behaviour appears to be unaffected by the inclusion of strat- 
ification and weak compressibility effects. This has been analyzed 
only for the case of constant soundspeed (i.e. constant disc back- 
ground temperature profile). In Section 4 we discuss these results 
and conjecture as to the disappearance of the SRI in this inviscid 
QHSG limit. 

The main results of this paper can be summarized here: (a) the 
SRI is also present in the LSB model equations, (b) we find that in 
the inviscid-QHSG limit only the doubly reflecting boundary con- 
ditions, i.e. no-normal flow at the two radial boundaries, seem to 
admit unstable normal modes, (c) the inclusion of stratification re- 
alistic for a disc (again, in the QHSG approximation) does not alter 
the SRI as obtained in previous investigations, and (d) there is an 
advectively conserved quantity in the BE which has the character 
of a potential vorticity; its existence implies both that this type of 
flow always has a continuous spectrum associated with it and that 
there exists certain conditions, including large Brunt- Vaisaila fre- 
quencies, in which the underlying flow behaves two dimensionally 
by advectively conserving the vertical component of its vorticity. 

Note that the terms channel flow and no normal flow are used 
interchangeably in this manuscript. 



1 EQUATIONS AND BOUNDARY CONDITIONS 

1.1 The Large Shearing Box and Boussinesq Equations 

There are a few formal asymptotic derivations for the set of equa- 
tions appropriate to the dynamics of a localized section of a rota- 
tionally supported disc (e.g. Goldreich & Lynden-Bell,1965). For 
this discussion we shall begin with the Large-Shearing Box (LSB) 
equations as they appear in Umurhan & Regev (2004), 



(d t - qCloxd v )p + V ■ (p b + p)u = 0, 
(dt — qQoxd y )u + u ■ X7u — 2flov' = — 



dxP 
Pb + p 



(1) 
(2) 



dyP 



(dt — q£loxd y )v' + u ■ Vv' + (2 — q)Vlou = ^—(3) 

pb + p 

(d t -qn xd y )w' + n'-Vw' = —^-^-, (4) 

Pb + p pb + p 

(d t - qn xd y )E + u ■ VE = (5) 
in which the total entropy is defined by 

Pb+P 



CVln 



(Pb + p)~ 



where Cy is the specific heat at constant volume and where 7 is the 
usual ratio of the specific heats at constant pressure and volumes. 



. These equations represent the dynamics taking place in a "large- 
shearing box" section near the midplane of a Keplerian disc rotating 
with radius Ro about the central star with the local rotation vector, 
f2(i?o)z. 

The above equations have already been non-dimensionalized. 
Time is scaled by the local rotation time of the box. All lengths 
have been scaled according to a length scale L which is compara- 
ble to the disc thickness. Pressures are scaled according to the local 
scale of the soundspeed, which is in turn based on some fiducial 
characteristic temperature scale. For further details see Umurhan & 
Regev (2004). The above equations, in which the vertical compo- 
nent of gravity is constant, are identical to the equations considered 
by Tevzadze et al. (2003). 

In the language of this paper, x corresponds to the shear- 
wise (radial) coordinate which is a small section located around 
the disc radius Ro, the azimuthal coordinate y is streamwise and 
z is the vertical coordinate corresponding to the normal direction 
of the original disc midplane. The velocity disturbances in the ra- 
dial, azimuthal and vertical directions expressed by the variables 
u' = {u , v' , w'}. These velocities represent deviations over the 
steady Keplerian flow. 

Slo, sometimes also referred to as the Coriolis parameter, is 1 
in these nondimensionalized units, meaning to say because time has 
been scaled according to the dimensional value of the rotation rate 
at Ro, i.e. Q(Ro), the local Coriolis parameter formally is equal to 
one. We retain this symbol in order to flag the Coriolis effects in 
this calculation. The local shear gradient is defined to be 



(6) 



in which Q(R) is the full disc rotation rate. For Keplerian discs the 
value of q is 3/2. The local Keplerian flow is represented here by a 
linear shear in the azimuthal direction, i.e. — q£loxy. 

The derivation of the LSB equations from the full-scale disc 
equations begins from the observation that the disc is in fact rota- 
tionally supported to leading order. This means that in steady state 
the disc radial force balance is, to leading order, between rotation 
and the radial component of gravity as emanating from the central 
compact object. The radial steady state pressure gradients provide 
corrections which are on the order of e 2 = H 2 /Ro, where H ia the 
characteristic height of the disk. Cold discs are those in which e is 
formally assumed to be a small quantity. Given the scaling and ex- 
pansion steps leading to the LSB this rotational support translates 
to saying that a radial buoyancy term in a rotationally supported 
disc is asymptotically smaller than the Coriolis effect by O (e 2 ) . It 
is for this reason that there is no buoyancy term in J2j. However, 
this rotational balance implies that the vertical structure of the disc 
is primarily characterized by the usual hydrostatic balance since, 
by the symmetries inherent to the geometry of the problem, there 
is no rotational support in that direction. It is for this reason there 
is an associated buoyancy term in J4j. It should be noted that in the 
context of the LSB, the steady state pressure and density functions 
are only functions of the 2 coordinate (see below). The vertically 
oriented gravitational field emanating from the primary compact 
object is linearly proportional to the distance from the midplane in 
the LSB system (again see next paragraph). This is a consequence 
of the asymptotic expansions leading to the LSB in which the cen- 
tral object's gravitational potential is expanded in a Taylor series 
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% Note that incompressibility has the mathematical effect of an infinite 
soundspeed. 



1 In Umurhan & Regev (2004), the entropy (heat) equation was written ex- 
plicitly in terms of the pressure and densities. There is, then, no substantive 
difference between these two expressions of the same conservation law 
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about the disc's symmetry plane (i.e. z = 0). Again, for details of 
this procedure see Umurhan & Regev (2004). 

The steady state density and pressure functions are given by 
Pb, pb- These relate to each other according to the aforementioned 
local hydrostatic equilibrium relationship, 

d z p b = -pbg(z), 

with g(z) — QqZ. The detailed solution for the steady states are 
then determined once something has been said relating the steady 
quantities. For simple theoretical investigations this comes in the 
form of a barotropic equation of state, namely, the statement that 
Pb = Pb(pb)- Note that in this paper, especially in Section 4, we 
explicitly assume that these steady quantities depend only on z and 
not on x. The dynamic pressures and densities, p and p, represent 
deviations about their corresponding steady state quantities. 

We want to gain some insight as to the effects that gravity and 
gradients of state quantities have on the local dynamics. We take an 
incremental approach towards this goal by considering a more sim- 
plified version of these equations. In more concrete terms, the linear 
theory i 1 15 i will show the presence of three "types" of temporal 
modes, (i) a pair of acoustic modes, (ii) a pair of inertial-gravity 
modes and (iii) an entropy mode (Tevzadze et al., 2003). Although 
the acoustic modes are interesting, we chose to consider the dynam- 
ics of a system in which the acoustic modes are effectively filtered 
out. 

To do this we invoke the Boussinesq approximation which, in 
this sense, we replace the continuity equation Q with the statement 
of incompressibility and we replace the entropy conservation equa- 
tion {5J with an evolution equation for the temperature fluctuation, 
8. All density fluctuations are set to zero accept the one associated 
with the buoyancy term in in which it is related to the tempera- 
ture fluctuation via 

p' = -a p 6; ^.-g^. 

In other words, a v is the coefficient of thermal expansion at con- 
stant pressure. This is the typical formulation of the Boussinesq 
approximation (Spiegel and Veronis, 1960). We furthermore posit 
that in this limit the basic state density profile is a constant and in 
order to distinguish this from a spatially varying density profile we 
denote the former with pb. The resulting model set of equations 
are similar to those assumed in the studies by Yavneh et al. (2001), 
Dubrulle, et al. (2005) (these being viscous studies) and Rudiger et 
al. (2005) (a cylindrical Taylor-Couette analysis). We have then, 
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(10) 




— qCloxdy 


)6' + u' 


• V0 + wd z T b = 0. 




(11) 



These are otherwise known as the Boussinesq equations (BE) in 
plane-Couette shear. The term TJ, is the way in which the basic 
state temperature profile varies with height in this model formu- 
lation of the disc system. Such a term would vary according to 
d z Tb = T zz z, in other words, the gradient of the basic state tem- 
perature has a linear dependence with respect to the disc height in 
which T zz is the parameter that controls the slope of this varia- 
tion. For situations in which T zz is negative, the atmosphere can be 



thought of as being classically buoyantly unstable which could lead 
to Rayleigh-Benard convection (see for instance, Cabot, 1996). 

The BE equations here are mathematically equivalent to the 
inviscid limit of the equations studied in D05. The only difference 
here is in interpretation. Whereas we follow a temperature pertur- 
bation, 8, and steady temperature profile, T b , they follow an entropy 
perturbation, denoted by h, and steady entropy distribution H. The 
two disturbance quantities are related to each other via 

ja p C v 



h = -- 



Pb 



(Also see Appendix [DJ. 7 is the ratio of specific heats and C v is 
the specific heat at constant volume. The equations 1711 li are the 
basis of the discussion in Section 3. The steady configuration of 
these equations which will be perturbed is 



0, 



p — p — constant. (12) 

The set I l!5i . including their corresponding analogous bound- 
ary conditions will be considered in Section 4. 

1.2 Boundary conditions 

There is no obvious choice of boundary conditions for this reduced 
set of inviscid flow equations. We consider all disturbances to be 
periodic in the y direction. Because these are equations meant to 
model what happens near the midplane of a circumstellar disc, we 
can consider periodic conditions in the vertical direction only under 
special circumstances (i.e. the constant Brunt- Vaisaila frequency 
approximation of the BE in Section 3, see below). In a more gen- 
eral sense we will distinguish, instead, between either varicose or 
sinuous modes. By varicose modes we mean to say disturbances 
which have even symmetry with respect to the z — plane in all 
disturbances except the vertical velocities, which have odd symme- 
try with respect to the 2 = plane. Sinuous modes have the reverse 
symmetry of the varicose modes. In situations in which a bound- 
ary condition needs to be specified on vertical boundaries of the 
atmosphere, we assume that there is no normal flow. 

The more troublesome of the boundary conditions has to do 
with what to say about the flow variables in the radial direction. 
This is because an injudicious choice of a boundary condition 
might incite disturbances of the fluid into instability by drawing en- 
ergy across the boundaries. It is our interest here to minimize this 
potential as much as possible. Of the myriad of possible choices, 
we see the following three sets of boundary conditions as ones that 
achieve this objective (and motivated further in Appendix^): (a) 
that the flow be confined between channel walls lying at x = 0,1, 
which means in practice that the radial velocities are set to zero 
there, i.e., no normal-flow conditions; (b) the flow has zero La- 
grangian pressure fluctuations (defined below) on both of the mov- 
ing radial boundaries and; (c) a mixture of these two conditions, 
for example, by requiring there to be no normal-flow on the inner 
boundary while there is no Lagrangian pressure fluctuation at the 
outer boundary. 

The vanishing of the Lagrangian pressure fluctuation on the 
undulating radial bounding surface S r in linear theory translates to 
requiring 

p + ixdxp = 0, 

where p' is the pressure fluctuation about the steady state. The po- 
sition of any particular radial surface, initially at rest at coordinate 
x, is denoted by £ x (y, z), and evolves according to its Lagrangian 
equation of motion (Drazin & Reid, 1984) 
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u'(x, y, z, t) = ^ = (d t - qCloxd v )£ x + v'd v £ x + wd z £, x .(l3) 

Because the steady state pressure configurations of both the BE 
(p) and the LSB (pb) are constant with respect to x, the condition 
simplifies to requiring 



P =0, 
atx = 0, 1. 



2 LINEAR DYNAMICS OF THE BOUSSINESQ 
EQUATIONS 



(14) 



In this inviscid limit there emerges a natural timescale defined as 
the Brunt-Vaisaila frequency, N. This time scale is defined through 
the product of the vertical temperature gradient and vertical gravity 
via 



N 2 = gj- b a p d z T b 



(15) 



which is, in general, a function of the vertical coordinate z. 
Throughout this study N is taken to be real (buoyantly stable). Lin- 
earization of 171 111 reduces to, 

(d t - qfl xdy) u - 2n Q v = -±-d x P, (16) 

l J b 

(dt-qCl xd v )v + n (2-q)u = ~j;d y P, (17) 



(dt — qSlaxdy) w 



(d t - qtioxdy) 9 = -N 2 w, 
d x u + d y v + d z w = 0, 



(18) 

(19) 
(20) 



where the temperature variable has been slightly redefined as O = 
gctpO/pb- pb is set to 1 from here on out. It is now to be under- 
stood that unprimed velocity expressions (i.e. u,v,w) represent lin- 
earized disturbances. 



2.1 A conserved quantity for linearized flow 

There exists a conserved quantity in these equations. Operating on 
< 1 61 by d y followed by operating I17i by d x and subtracting the 
result reveals 



(dt — qiloxdy) (d x v — d v u) = Q (2 — q)d z w, 



(21) 



where the incompressibility condition has been used. The term on 
the LHS of this expression is the vertical vorticity, i.e. £ = d x v — 
d y u. With a similar tack one can multiply {T5J by Q (2 - q)/N 2 
and then operate on the result with d z to get 



. d go(2-g) n 
(d t - qfloxdy) [ — — 9 



-U (2-q)d z w. (22) 



Adding the results together yields a general conserved quantity of 
linearized flow of this type: 

(d t - qtl xd x ) E = 0, (23) 
where 

The quantity S can be thought of as a generalized potential vor- 
ticity for this type of flow whose analogous quantity is discussed 
in Tevzadze et al. (2003). The conservation of S immediately im- 
plies that there always exists a continuous spectrum (see Schmid 
& Henningson, 2001, Tevzadze et al., 2003) for this type of phys- 
ical system. Note also that the system of linearized equations (i.e. 



1161201 is third order in time. One may suppose that there are three 
independent normal modes for any given set of quantum numbers 
of the system (see below), however, given that there exists a con- 
served quantity, together with its associated continuous spectrum, 
it means that there are at most only two normal modes for any given 
quantum number set. 

Inspection of S shows that disturbances behave in a quasi two- 
dimensional fashion in some limits. One of these is when q = 2, 
that is at the critical Rayleigh condition (Drazin & Reid, 1984): 
it follows that the vertical vorticity is conserved by the flow. The 
second of these is to notice that if the temperature fluctuation re- 
mains an order one quantity as N 2 gets large then the flow again 
exhibits quasi two-dimensionality with the vertical vorticity be- 
ing conserved. We reflect upon the consequences of this conserved 
quantity some more in the Discussion. 



2.2 Constant N 2 

Since part of the purpose of this work is to further develop some 
amount of intuition about the dynamics of such disc environments 
primarily through analytical means, it will be more tractable for us 
to first treat the vertical gravity component g(z) to be 



g(z) = gosgn(z). 



(25) 



The non-dimensional constant go is technically arbitrary. In a sim- 
ilar vein we approximate the steady state temperature gradient by 
saying 



d z T b = T z sgn(z) 



(26) 



in which T z is another non-dimensional parameter. The conse- 
quence of this is that TV 2 is a constant for z 7^ 0, and is zero at 
2 = 0. At this stage, these assumptions are qualitatively no dif- 
ferent than what has be done in Yavneh et al. (2001), D05 and 
Shalybkov & Rudiger (2005), although we take a more realistic 
interpretation of a constant N 2 (and see below). From here on out 
we set pb = 1. Additionally, we restrict analysis of the dynamics 
to z > and keeping in mind that modes are considered to have 
either sinuous or varicose spatial character in the vertical. 
We write general solutions into the form 




cos f3z 
sin f3z 



while for the other variables 



e. 



— sin /3z 
cos fiz 



e + c.c, 



e + c.c, 



(27) 



(28) 



The terms above in the curly brackets represent varicose distur- 
bances while the terms below are the sinuous disturbances. In this 
sense the "quantum numbers" of the system are given by a, j3 (vari- 
cose or sinuous) and a radial overtone number (if there are more 
than one) subject to solution of the normal mode boundary value 
problem below. 2 . Note that because this is a single Fourier expan- 
sion we restrict our considerations to < a < 00 together with 



2 The use of quantum numbers should be considered only in terms of con- 
ventional nomenclature. There is no real quantization in the horizontal and 
vertical directions per se since we allow these quantities to take on any value 
from the continuum of real numbers 
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< P < oo. Insertion of 1271281 into the governing linear equa- 
tions gives, 



i (uj — qfloxa) u af} — 2f2ow Q(3 = —d x P a 

i (ui — qQ a xa) v a/3 + tto(2 — q)u al3 = —iaP a 

i (lu — qfl xa) w a „ = —/3P a/3 

i (uj — qVloxa) Q/3 = —N 2 w < 

P W a = °. 



e 



(29) 
(30) 

•p. (3D 
(32) 
(33) 



It turns out that it is much more tractable to consider the lin- 
earized normal-mode behavior in terms of equations describing the 
pressure fluctuation and radial velocities. This is entirely analogous 
to what was done in D05 and the following equations should be 
compared to the ones quoted in D05 as Equations (21-22). 



f3 2 o 2 



N 2 



I 2 

( CT 



iP 



2\ 



T(). r 



= (ad x iP a 



+ Qq(2 — q)au a 



2£l aiP 



(34) 



where, for the sake of compact notation we use the expression 
a = lj — qfloax. The epicyclic frequency to 2 is equivalent to the 
expression 2(2 — q)Q 2 - 

We are mainly interested in analytically expressible solutions 
to the above set of equations. To achieve this in an asymptotically 
rigorous manner the following scalings seem natural: when the hor- 
izontal wavenumber is small, it follows that the frequency scales 
similarly. Using e to measure this smallness it follows, 



eax, 



eoJi + ■ 



To lowest order it also follows that a = eci + 
and velocities are consequently expanded by 

P„ B = Pa + 2 P 2 + ■ ■ ■ 



The pressure 



u a/3 = eui + e U3 H 

Implementing these expansions into the governing equations J34I 
yields at lowest order in e a single equation for the pressure pertur- 
bation, 

(wi - qQoaix) (d 2 x - F 2 P 2 ) P = 0. (35) 
Normal-mode type solutions to 1351 are, 

Pq = A cosh k F x + B sinh k F x, (36) 

where the Froude-wavenumber, k F , is defined as k F = 
uj 2 f3 2 /N 2 — p 2 F 2 . The epicyclic-Froude number is denoted by 
F e , This mathematical structure of I35i is identical to the operator 
describing the evolution of plane-Couette disturbances in a channel 
(e.g. Case, 1960). 

We compare the analytical solutions generated here with nu- 
merical solutions generated for the boundary value problem defined 
by the un-approximated full linearized equation set I34i . A sec- 
ond order correct (in the x direction derivatives) Newton-Raphson- 
Kantorovitch (NRK) relaxation scheme on a grid of approximately 
1000 to 2000 points is used for the verification. Relative conver- 
gence was checked by doubling the size of the domain. Eigenvalues 
are determined with errors that were no more than O (5 x 10 -6 ). 
As such the eigenvalues generated asymptotically are considered to 
be valid in all cases where normal modes exist. 

Because this is a relaxation method reasonably good initial 
guesses are required, both in the eigenfunction and the eigenvalue, 
in order to accurately obtain an answer. When the initial guesses 
were far off from the actual solution the scheme admitted solu- 
tions belonging to the continuous spectrum. This occurs for all 



sets of boundary condition but is the only solution possible for the 
case of fixed pressure conditions (see below). Therefore we discuss 
the general features of the continuous spectrum in Section 12.2.21 
Nonetheless, we find that it helps to avoid the continuous spectrum 
if the initial eigenvalue guess is set so that Im(w) 7^ 0. The pos- 
sibility of the solution jumping onto a random continuous mode 
solution is satisfactorily bypassed in this way (also see below). 

In what follows we consider the discussion for each of the 
three boundary conditions. We finally note that the resulting math- 
ematical structure resulting from this asymptotic limit is closely 
similar to the WKB analysis done in D05. Whereas in this study 
the small parameter is the horizontal wavenumber, in D05 the small 
parameter is the shear term qQo (denoted by 'S' in their equation 
3). In this sense, their asymptotic form is valid for small values of 
the shear while ours is valid for small values of a. 



2.2.1 No normal-flow conditions 

The wall conditions, i.e. that u a/3 
terms of the pressure the requirement 



at a: = 0, 1, becomes in 



-iioAPo + 2n QiPo) = 0, at x = 0, 1, (37) 



at lowest order. Implementing these conditions and a little algebra 
gives a dispersion relation for uj\ 



cngOo I - ± — - — A]/ 2 
' 1 2 2k F q 



A F = 16 + hi 



8k f q 
tanhfc^ 



(38) 



As the dispersion relation clearly indicates, if A F < then there 
appears a pair of complex modes, one which grows and one which 
decays. When A F > there are two propagating modes oscillating 
with no overall growth in amplitude. The character of the stability 
is dictated only by the two parameters q and k F . The limit where 
A^ 2 — > (i.e. k F — > 00) reveals that uj\ — > 0, qaiQo- 

The striking feature of this general solution is that there ex- 
ists a band of vertical wavenumbers for which a stable/unstable 
solution exists. In Figure Q we plot this dispersion curve for the 
case q = 3/2. The plot shows a band in Froude-wavenumber 
within which the stable/unstable pair exists. Recall that the Froude- 
wavenumber is really the vertical wavenumber scaled by F e . The 
bifurcation into the stable/unstable pair occurs when the frequen- 
cies of the two modes become the same. In the case depicted in 
the figure (q = 3/2), when the frequencies merge the instability 
emerges near kF ~ 2.1. Until about kF ~ 3, where the instability 
vanishes, the frequencies remain the same. 

The boundaries of the unstable band for general values of q 
may be inferred from the expression for the growth rate: this means 
determining the function q±(kp) that satisfies A/(g, kp) = as 
defined in <38> . The two functions are g_ = 4tanh(kF /2) /kF and 
q+ = 4coth(kp/2)/kF. We see that the band structure for the 
instability range persists until q = 2, which happens to also corre- 
spond to the Rayleigh instability line. Beyond q = 2 the instability 
range is bounded from below by zero vertical wavenumber but it is 
still bounded from above by a finite j3. As the shear becomes weak, 
the band of unstable modes gets correspondingly thinner. Note that 
this instability disappears in the non-shearing limit, that is when 
q — > 0. Figure|2|graphically summarizes these results. 
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2.2.2 Fixed pressure conditions 

The boundary conditions on the lowest order pressure conditions 
becomes 

Po = 0, (39) 

at both boundaries x = 0,1. Given the form of the underlying 
equations it turns out that there is no normal mode solution possi- 
ble. The numerical procedure admits solutions, however, these are 
always modes associated with the continuous spectrum of the lin- 
earized system. They are not true normal modes in the usual sense 
because they exhibit a discontinuity in some quantity: here being in 
the horizontal velocity, v, and manifesting explicitly as a step in the 
quantity d x P. Discontinuities of this sort, referred to sometimes as 
singular eigenfunctions, are typical features of modes associated 
with a continuous spectrum (Case, 1960, Balmforth & Morrison, 
1999). In all cases, the location of the discontinuity is at some value 
of x = x c which is the location of the critical layer, in other words, 
the place where the quantity lui — qQ,oa\x is zero. This means (and 
this is verified numerically) for given values of a,q and Qq there 
will be a continuum frequencies, ui c (a, q, f2o), existing between 
and qQ.oa. 

Figure lBTl displavs examples of this continuum mode together 
with an analytic representation of the continuum mode <B3> devel- 
oped in AppendixlBl 



Frequency Response for q = 3/2 




vertical wavenumber: kp = F g pS 

Figure 1. Asymptotic dispersion relation, £Ji. q = 3/2 for channel flow 
disturbances. The dashed curve is the growth/decay while the solid line is 
the frequency. The bifurcation into the stable/unstable state begins and ends 
near kp ~ 2.1, 3 respectively. In this range the frequency of both modes 
are identical. 



Stability Boundaries 



2.2.3 Mixed conditions 

We define terms by identifying mixed-A boundary conditions with 
zero radial velocities at x — and zero pressure perturbation at 
x = 1, while mixed-B boundary conditions indicate zero pressure 
fluctuations at x = with zero radial velocity perturbations at x — 
1. Both boundary conditions yield single normal mode solutions 
Consequently for mixed-A conditions the frequency response is 



lui — lloqai tannfc^, 

while for mixed-B conditions we have 

2^0 Qfi 



Wl 



■ tanh fc„ 



(40) 



(41) 



In Figure|5|we plot sample eigenfunctions for all boundary condi- 
tions we considered in these sections, as well as comparisons be- 
tween the analytic and numerical solutions obtained. 

It is important to mention that when discrete normal modes 
have frequencies which sit in the continuous sea, i.e. when < 
uj < qQ,oa it becomes challenging for the numerical method to 
not mistake it with a mode belonging to the continuous spectra. To 
circumvent this possible ambiguity (and to properly numerically 
verify this limit) we follow modes in j3 starting initially with values 
of the discrete normal mode frequency which is beyond the con- 
tinuum sea. In this way, discrete normal mode solutions are easily 
found and, using these as a starting point, one may incrementally 
move into the realm where normal modes exist within the contin- 
uous sea. This is depicted in Figure |4| In all boundary condition 
cases investigated, we successfully trace the discrete mode spec- 
trum. 



2.3 The QHSG approximation: vertically varying A?" 2 

The asymptotic results of the last section hints toward a tractable 
approach in evaluating the generality of the SRI under a variety 
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\ Rayleigh Stable i 




q + = 4coth(k p /2)/k F 


Stable ! 


q = 4tmhfly2)/k^ " — 
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0.5 1 1.5 2 2.5 3_ J.5 4 4.5 5 5.5 6 
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Figure 2. Plotted is the stability boundary for channel flow disturbances in 
the limit where a is small. The shaded region corresponds to the instabil- 
ity regime. The bottom portion of the instability region is bounded by the 
curve g_ = 4tanh(.F e /3/2)/.F e /3 while the upper portion is bounded by 
= 4coth(F e /3/2)/_F e /3. The marginal stability values for Keplerian 
flow, i.e. q = 3/2, are labeled by dashed lines. The Rayleigh unstable line 
is indicated, i.e. q > 2. 



of conditions. The limit where the horizontal wavenumbers (i.e. 
a) are small suggests that there exists well-posed reduction of the 
governing equations of motion 181 1 1 i . With the azimuthal scales 
of disturbances scaling as 0(1/ a) it followed that the temporal 
disturbances scale as O (a) and, as such, it implied that the radial 
and vertical velocities similarly scale as O (a) while the pressure, 
the temperature fluctuation and the azimuthal velocities scale as 
O (1). We therefore propose that when the horizontal scales are 
large compared to the corresponding vertical and radial ones the 
following scalings hold, in general, with respect to quantities and 
operators of the system: 



dt,d y ,u,w 



O(a) 



d x ,d z ,p,v, 



0(1). 



(42) 
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a = 0.35, p = 1 
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a = 0.1,0 = 0.35 



Normal Mode Dispersion Relations 
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- ' Channel Walls 
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0, 
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pb 



d x u + d y v + d z w 
+ (a 2 ) 

(d t — q0.oxdy)v + u ■ Vi> 

+ (a 2 ) 
{d t - qn xd y )6' + u • V0 = -wd z T b 



Pb 
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Figure 3. Assorted eigenfunctions for normal mode solutions are plotted for 
N = 0.4. In the first column shows the solutions for the parameters a = 
0.35, /3 = 1 while the second column depicts solutions for a = 0.1, /3 = 
0.35. Each graph shows the behaviour of the indicated quantity for the three 
sets of boundary conditions: channel walls (dashed-dot), Mixed-A (dashed) 
and Mixed-B (solid). Because for these parameter values the resulting u> are 
real, the real parts of the pressures and horizontal velocities are zero while 
the imaginary parts of the radial velocities are zero. 



This means that in the limit where the azimuthal scales are long, we 
have the following effective reduction of the nonlinear equations of 
motion flTTV 



(43) 
(44) 

(45) 

(46) 
(47) 



The above set is similar to the quasi-geostrophic, quasi-hydrostatic 
approximation used in the study of atmospheric flows (e.g. Ped- 
losky, 1987, Salmon, 2002). Whereas in the terrestrial analog full 
quasi-geostrophy involves retaining the inertial terms in the radial 
momentum equation, here they are absent (cf. !44> . and it is for this 
reason we consider the above set of equations to be a sort of semi- 
geostrophic limit. The semi-geostrophic nature of this set shares 
some similarities with the so-called elongated-vortices equations 
derived in Barranco et al. (2000). The set presented here differs 
from that work in that the elongated- vortex equations do not make 
the hydrostatic approximation as is a natural and necessary conse- 
quence here. 




1 1.5 2 2.5 3 3.5 4 
Froude-wavenumber: kj_ = F g p 

Figure 4. Normal mode dispersion plots for the three sets of boundary con- 
ditions considered. The frequency is shown scaled by qQ,$ot\. Plotted are 
the asymptotically determined dispersion relationships. The shaded region 
represents the continuous sea, i.e. modes associated with the continuous 
spectrum. Every point in the shaded region has a continuous spectrum mode 
in it. The strategy for following the discrete modes into the continuous sea 
using the numerical method begins with locating a discrete mode outside of 
the sea and then following the dispersion branch into the sea (as indicated 
on the plot with arrows). 



The power in this reduced set of equations, aside from exactly 
reproducing the asymptotic limit explored in the previous section, 
is that it allows one to investigate the effects that a position de- 
pendent function of gravity and background state temperature gra- 
dient, i.e. g(z) and d z T b , have on the SRI. In this sense, unlike 
the approximation utilized in Section |2~2| we relax the condition 
that g and d z T b are constants and let them be general functions 
of 2. It therefore means that the Brunt- Vaisaila frequency is now 
z-dependent. 

When specific forms are considered here we assume that these 
quantities are simply proportional to z, that is to say, 



g(z) = Q z; 



d z T b = T z: 



(48) 



The constant T zz sets the severity of the background temperature 
gradient (see Section 1.1). 

We linearize the set 1431471 about the quiet state u = v — 
w = = 0, p = constant. Disturbances are denoted with primes. 
A little algebra shows that these equations may be simplified into a 
single one for the pressure perturbation: 



(dt — Qqxdy) 







dp' 



dzN 2 (z) dz 



+ 



a p 

dx 2 



= 0, 



(49) 



where the z-dependent Brunt- Vaisaila frequency is given as 

N 2 (z) = j- b a p gd z T b = j-nlf zz a p z 2 . 

This all also means that we can consider a z-dependent Froude 
number according to 

, .2 1 



N 2 (z) 



F 



(50) 



because all variables and quantities have been non- 
dimensionalized, the Froude-number scale F e should be con- 
sidered in parallel to the Froude-number F e treated in Section 
2.2. 
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The expression inside the square brackets of 1491 looks analo- 
gous to the potential-vorticity of atmospheric flows. When JV 2 (z) 
is a constant then the expression within the square brackets exactly 
recovers the asymptotically valid governing equation for the pres- 
sure perturbation in <35> . 

Separable solutions are assumed of the form 



V 



U(z)P (x)e 



iuit-\-ia.y 



+ C.C. 



(51) 



and applied to the governing equation 1491 . This results in two 
ODE's, one for the vertical structure function and one for the ra- 
dial structure function, 



dz e dz 

d 2 p. 



dx 2 



-kin, 



klPo 



(52) 
(53) 



The separation constant for this procedure is k F . We notice imme- 
diately that the equation for the radial structure function is math- 
ematically the equivalent to 1351 . Since the boundary conditions 
and associated relationships are identical in this asymptotic limit, 
it follows then the same stability properties that was determined 
in Section l2~2l carry over here to this particular example of a z- 
dependent function of N 2 if the allowed values of k F are real. 
In remaining consistent with the terminology introduced in D05, 
modes for which k F are real are referred hereafter as exponential, 
or as e-modes, since this describes the quality of the radial structure 
function that results from solving 1531 . By contrast, modes in which 
the k F are imaginary are referred to as oscillating or as o-modes, 
(again see D05) and these in principle will have different stability 
properties than those determined for the e-modes. 

Thus the task that remains is to determine the allowed values 
of the separation constant. The analysis below indicates that that the 
only types of disturbances permitted for TV 2 (z) ~ z 2 are e-modes 
on a finite domain. 



2.3.1 Exponential modes 

Though these may be artificial, for the sake of simplicity and com- 
parison we consider only the vanishing of the normal velocities at 
the boundaries z — ±1 (e.g. Barranco & Marcus, 2005). When as- 
suming the specific forms for g and d z Tb as in 1481 we find two 
possible solution forms to 1521 



n(*) 



(varicose) 
(sinuous) 



J -i V 2 ^ 



(54) 



where the symbol J denotes the Bessel function of the first kind. 
The expression k F /F e can be considered to be parallel to the verti- 
cal wavenumber /3 introduced and treated in Section 2.2. A Taylor 
Series expansion of these solutions near z = verifies the even 
(odd) symmetry of the varicose (sinuous) solutions. Given these in- 
herent symmetries we are left with the task of setting to zero w' at 
z — 1 which, given the relationships between 146147 K is equivalent 
to setting d z p' — there. Given the general properties of Bessel 
Functions (Abramowitz & Stegun, 1972) the set of k F values that 
satisfy the boundary conditions are always real. In Figure|5|we dis- 
play the functions n(z) for the first three values of the separation 
constant k F . It should be noted that the asymptotic expansion of the 
solution forms presented above are characterized by an amplitude 
function which grows as <Jz for \z\ — > oo. 
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Vertical Eigenfunctions: First 3 modes 
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Figure 5. The vertical eigenfunctions Tl(z) for e-modes are depicted for the 
first three modes (varicose and sinuous) that satisfy the boundary condition 
that the vertical velocity is zero at z = ±1. Due to the inherent symmetries 
in the modes, the behaviour of the function is depicted for z ^ 0. 



2.3.2 Oscillating modes 

Consideration of o-modes starts by redefining the separation con- 
stant k F according to 

k F itz F , 

where k f is real. Therefore, as in the previous section there are two 
independent solutions to I52i with k F so defined given by 



1 K p 

4 V2 F e 



z*I 



1 K F 2 
2F/ 



respectively corresponding to sinuous and varicose disturbances. 
However, Modified Bessel Functions as these do not have zeros 
for real values of the ratio K F /F e . It means that it is not possible 
to satisfy the same sort of boundary conditions considered in the 
previous section (e.g. w = at z — ±1) which, in turn, means that 
o-modes are not allowed solutions on a finite vertical domain with 
no-vertical flow boundary conditions. 

However, a consideration of the problem of J52I on an infinite 
domain (i.e. z — > ±oo) where, instead, we require boundedness 
and asymptotic decay of all quantities as \z\ — + oo, suggests that 
the solution for H(z) could be the following, 



'"--» = ■ 1 ' v - ( \y* 2 



(55) 



where K„ (x) is the Modified Bessel Function of the Second Kind 
(Abramowitz & Stegun, 1972). Indeed an asymptotic expansion of 
the leading order behavior of /C 3/4 (z 2 ) in the large argument limit 
shows that it behaves like a Gaussian, i.e. like exp — |^z 2 , for 

large z 2 . On its surface this behaviour seems to satisfy the require- 
ments on the functions as z — » ±oo. Although the above conclu- 
sion is correct for z — > oo, to extend this conclusion as z — > — oo 
based on the above representation would be wrong. In fact, given 
the functional form <55t . it becomes a matter of subtlety as to how 
one must cross the point z = 0. In Appendix Iclit is shown that the 
behaviour of 1551 for z < is 
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(- 



/C, 



2 F r ^ 



- 2i i {m- 



(56) 



Inspection reveals exponential divergence as z — > — oo since T„ (a;) 
behaves exponentially as i — > oo. It appears that this analysis 
indicates that there are no bounded solutions possible for 11(2;) 
and, consequently, it implies there are no o-modes permitted when 
TV 2 (2) ~ z 2 in the context of the BE model. 



3 LARGE SHEARING BOX EQUATIONS: THE QHSG 
APPROXIMATION AND LINEARIZED DYNAMICS 

Given the clues revealed by using the QHSG for the BE, we con- 
sider the same scalings expressed in 1421 and apply them to the full 
LSB equations (1-5). The major departure here, of course, is that 
disturbances are now not incompressible. One scaling relationship 
is to say that the density and pressure variables are of comparable 
scale, i.e. O (p) = O (p) ~ 1. Therefore the QHSG reduction of 
the nonlinear LSB becomes 



(9, 

= 2fi « 



qQ,oxd y )p + V ■ (pb 
d x p 



p)u' = 0, 



p b + p 



+ 0(a 2 ) 



{dt — qQ.oxd y )v' + u ■ Vit + (2 — g)OoU = — 



d y p 
Pb + P 



= -a zP - pg{z) + O (a 2 ) , 
{d t - qQ xdy)E + u ■ V(E 6 + E) 



0. 



(57) 
(58) 

(59) 

(60) 
(61) 



where we have introduced the basic state entropy Ef, and its dy- 
namically varying counterpart E which are defined by 



E b = ln^, 
Pb 



In 



1 + -*- 

Pb 

(' + *)' 



(62) 



where 7 is the the usual thermodynamic ratio of specific heats. Lin- 
earizing I57I6H and sorting through the algebra (see Appendix ID! 
leaves us with a single master equation for the pressure perturbation 



{dt—qftoxdy) 



^d 2 p + 9 *J^ (Jp + d -P) + d *P 



The generalized Brunt- Vaisaila frequency is defined by 



N 2 



9 r, , Pb 

—a z In -=■, 

7 Pi 



0.(63) 



(64) 



(65) 



while the (nondimensional) adiabatic soundspeed c is defined by 

Pb 

<63> is the LSB equivalent, in this QHSG limit, of a potential vor- 
ticity for a local section of a circumstellar disc. Comparing this 
equation for the potential- vorticity with the analogous one for the 
BE in 1491 reveals some differences between them being, namely, 

uj 2 dp d uj 2 
g dz 



-i-p. 

dzN 2 c 2P 



The first of these is associated with the time rate of change of the 
density fluctuation in the continuity equation. This is explicitly ab- 
sent in the BE due to the assumption of incompressibility. The sec- 
ond of these is associated with the generalized entropy fluctuation 
and is inversely proportional to the soundspeed. This term is ab- 
sent in the Boussinesq Equations because the assumption of incom- 
pressibility is equivalent to the interpretation that the soundspeed is 
infinite. 



Having N 2 < is equivalent to the Schwarzschild condition 
for buoyant instability (Tassoul, 2000). As before, we assume that 
the atmosphere is stable to buoyant oscillations (iV 2 > 0). How- 
ever we must also say something about the soundspeed c: for the 
sake of this discussion we will assume that it is a constant with 
respect to z, that is, we assume the atmosphere is isothermal. We 
proceed toward determining normal mode solutions of the expres- 
sion inside the square brackets of I63i . 



lo 2 dp d 
g dz dz iV2 



' (9 ^ d P\ 



cPp 
dx 2 



: 



(66) 



Assuming separable solutions of the form J5 1 1 we find, once again, 
the following two problems to solve: 



g dz dz N 2 \c 2 dz J F 

d 2 p 



dx 2 



= KPo 



(67) 



(68) 



The separation constant k F is the same as before. 

Because the equation for Po is the same as in the BE model, cf. 
<53t . it immediately follows that the same stability properties that 
was determined for the BE apply here too if the set of k F values 
are all real (i.e. e-modes). For this study we restrict our attention to 
finite vertical domains (see below). It means, then, that the task that 
remains is to determine the eigenvalues of k F by seeking solutions 
of 1671 subject to the boundary condition that the vertical velocity 
vanishes at z — ±1. In the LSB, this condition amounts to setting 



N 2 \c 2 



(69) 



at 2 = ±1. (see Appendix D). Also, as before, we consider solu- 
tions to n that are either sinuous or varicose. 

We emphasize that we are restricting our attention here to so- 
lutions on a finite 2 domain. This is because attention needed to 
treat such problems on an infinite domain is challenging and it is, 
thus, outside the scope of this current work (see Discussion). 

Aside from very special values of the parameters, there are no 
simple or analytically tractable solutions to the ordinary differential 
equation posed by I67i 3 . Therefore we numerically solve for this 
equation and k F using a fourth-order variant of the NRK scheme 
discussed in the previous section. We use a grid of 300 points which 
lets us determine solutions up to machine accuracy (i.e. an error 
of less than 10~ n ). The solutions were all normalized by setting 
n = 1 at 2 = 1. We verify the robustness of the numerical scheme 
by using it to solve the simpler equation I49i and comparing the 
numerically generated results against the exact solutions 1541 . 

There are two parameters that govern the solutions. The first 
of these is the scale measure, F E , of the height-dependent Froude- 
number 

, „ 2 
Ft 



N 2 



Given that this atmosphere is isothermal this Froude-number scale 
is measured by the parameter 



u 2 y 



ng(7- 



For a medium dominated by molecular hydrogen 7 w 7/5. It 
means that in a Keplerian flow F 2 w 7/2. The second parameter is 

3 General solutions of this equation are linear combinations of hypergeo- 
metric functions which require numerical evaluation anyway. 
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Figure 6. Vertical wavenumbers and eigenfunctions for the fundamental 
mode solutions of 1671 for selected values of H and F 2 . 



the relative measure of the vertical scale height of the atmosphere 
defined by H and given to be 



H 2 



( >- 



in which 7£ M is the non-dimensionalized gas constant for the given 
composition, T is the non-dimensionalized temperature of the at- 
mosphere. This quantity is essentially the same as the classic e- 
parameter governing thin-disc theory (Shakura and Sunyaev, 1973, 
Lynden-Bell and Pringle, 1974) When H is small, the atmosphere 
is very shallow and, consequently, cold. 

The solutions that we scan all show that the k F values are al- 
ways real indicating these are e-modes (cf. Section 2.3). We were 
unable to find purely imaginary solutions for k F on this finite do- 
main. Thus it implies that the stability properties determined for 
e-modes (i.e. Section 2) carry over to here too and that this system 
does not support o-modes on this finite-domain. 



4 SUMMARY AND DISCUSSION 

4.1 The QHSG and the persistence of the SRI with height 
dependence in N 2 

We have achieved here an extension of the SRI to models which 
takes into consideration the vertical structure of the physical envi- 
ronment. One of the departures taken here from previous work is to 
explicitly include the effects of a vertically varying Brunt- Vaisaila 
frequency. The results of the previous sections shows that the SRI, 
under channel-wall boundary conditions, persists unaltered irre- 
spective of the model equations considered (i.e. either the BE or 



the LSB) or type of mode (i.e. e-mode or o-mode) so long as one is 
in the inviscid-QHSG asymptotic limit. This is not to say that if one 
relaxes the restriction of long horizontal length scale disturbances 
(i.e. small a) then this instability will not continue - we merely 
mean to say that its existence under those conditions remains open. 
It does seem likely, however, given the pattern of the presence of 
the SRI in D05, that it will do so also in the O (a) ~ 1 case too. 

The advantage of the inviscid-QHSG approximation em- 
ployed here is that an analytical analysis of normal-modes is possi- 
ble through a separation of variables. In general cases where both g 
and the other state variables like pt and pt are functions of the ver- 
tical coordinate z, the resulting linear equations and mode structure 
are generally non-separable. This makes for assessing the normal- 
mode behaviour of disturbances to be challenging (at best) although 
not impossible, as has been shown here. In this sense it appears 
that this is a reasonable peek into situations where complicating 
background structure in both the radial direction (here the shear) 
and vertical direction (here gravity and other state variables) can be 
taken into account together. 

It is also worth noting that in the inviscid QHSG limit of both 
the BE and LSB the radial structure of the eigenmodes (be they ei- 
ther exponential or oscillatory) are unaffected by the vertical strat- 
ification: in other words, the radial eigenfunctions are always the 
same. As a result, the stability criterion turns out to be insensitive 
to vertical stratification in the state variables of the system. 

The QHSG limiting process is achieved by looking at distur- 
bances with horizontal length scales that are large compared to the 
other dimensions which, therefore, implies that the associated time 
scale of disturbances can scale with proportional smallness. From 
the framework of the the LSB equations, the implication is that one 
recovers the part of the dynamics of inertio-gravity waves modified, 
to some extent, by the effect of weak compressibility (dilatation) 
and a finite sound speed. The inclusion of these effects is not to say 
that some facet of acoustic disturbances are recovered in this limit: 
this is precisely because the time scales associated with acoustics 
are much shorter than the time scales explored here. In this sense 
this asymptotic limit naturally filters out direct acoustic effects and 
preserves the dynamics of disturbances that would usually be as- 
sociated with Rossby waves in geophysical flows (Pedlosky, 1987). 
This is also not meant to imply that acoustic effects are unimportant 
(an issue that is far from settled), it merely means that this limit is 
effective at isolating the dynamics of these waves. 

4.2 On the nature of e-modes and o-modes with height 
dependent N 2 

D05 showed that there are two types of SRI modes which are char- 
acterized by the quality of the mode's radial structure function: 
either exponential or oscillatory. Because both the vertical back- 
ground temperature gradient and component of gravity are constant 
D05 showed that both e-modes and o-modes are supported for such 
a model atmosphere. 

For example, e-modes, which have vertical structure functions 
which are sinusoidal with respect to the vertical coordinate, are per- 
mitted in D05 even if the atmosphere extends mathematically to in- 
finity. Because N 2 is constant with respect to z, the corresponding 
vertical structure functions have no envelope growth or decay with 
respect to the vertical coordinate. As such, such solutions satisfy 
reasonable expectations of boundedness as the atmosphere extends 
indefinitely. 

By contrast the analysis of the QHSG limit of the BE equa- 
tions, in which N 2 depends on z quadratically, shows that the re- 
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suiting vertical structure functions for e-modes to have envelope 
structures which grow (~ \Zjzj). This fact makes it impossible to 
meaningfully impose boundary conditions or boundedness condi- 
tions on solutions as \z\ — > oo. The analysis performed here seems 
to indicate that e-modes are restricted to BE systems involving fi- 
nite vertical domains with height dependent iV 2 . 

O-modes in model atmospheres with quadratic dependence of 
N 2 with respect to z are not allowed. In atmospheres with con- 
stant iV 2 o-modes are admitted on account of the exponential de- 
cay of the vertical structure function. For N 2 (z) ~ z 2 (cf. l2,3,2l 
it appears there is no way to construct bounded solutions in the di- 
rections z — > ±00 simultaneously. It was also demonstrated that 
o-modes are ruled out on finite domains. 

The conclusions regarding the existence of o-modes is mainly 
based on the analysis of the BE with N 2 (z) ~ z 2 in the QHSG 
limit. We showed also in Section 3 that similar conclusions seem to 
hold for SRI e-modes and o-modes in the isothermal-QHSG limit 
of the LSB model set when considered only on a finite domain. It is 
not entirely clear how the existence properties of SRI modes are af- 
fected when: (i) the domain mathematically extends to infinity, (ii) 
the soundspeed varies with height. These questions are reasonable 
points of departure for further investigation. 



4.3 The absence of the SRI for non-reflecting boundaries 

The troubling aspect of this investigation is that when one consid- 
ers boundary conditions other than no-flow conditions on the ra- 
dial boundaries, the instability appears to vanish in the asymptotic 
limit considered. When the Lagrangian pressure is zero on both 
radial boundaries the analysis predicts that there are no normal- 
modes with the time scales assumed and, furthermore, there are 
only modes associated with the continuous spectrum. It is probably 
safe to conclude that for this type of boundary condition, that there 
are no normal-modes whose frequencies have magnitudes on the 
order of or greater than the small scaling parameter (a, the hori- 
zontal wavenumber) of the limit explored here. In situations where 
there is a mixture of no-flow conditions on one boundary and no 
Lagrangian perturbations on the other, only one normal-mode is 
admitted by the system which propagates with no growth or decay. 

The circumstance encountered here shares a number of simi- 
larities with the Eady problem of baroclinic instability in geophys- 
ical shear flows (Eady 1949, Criminale & Drazin, 1990). Drazin 
and Reid (1994) show that the Eady problem is essentially equiv- 
alent to the stability of inviscid plane Couette flow (pCf) subject 
to boundary conditions where the pressure perturbations are fixed 
on the channel walls, instead of the usual condition in which the 
normal velocities are set to zero. As Case (1960) showed, invis- 
cid pCf flows with no-normal flow boundary conditions admit only 
continuous spectrum modes and no discrete modes. By contrast, 
Criminale & Drazin (1990) showed that the Eady problem has, in 
addition to the continuous spectrum, a number of discrete modes 
present (which are possibly unstable under suitable conditions of 
the disturbances) when disturbances in inviscid pCf flow have fixed 
pressures at the boundaries (Criminale et al. 2003). 

For the SRI investigated here, an entirely analogous situation 
occurs: the linear operator governing the system is mathematically 
equivalent to the one characterizing inviscid pCf, but, the variables 
and stability characteristics are interchanged. Whereas in inviscid 
pCf the operator operates on the radial velocity (wall-normal) in 
the SRI case here it operates on the pressure perturbation. Thus the 
inviscid pCf problem admits normal modes (no normal modes) for 



fixed pressure (no normal flow) boundary conditions while the SRI 
problem admits normal modes (no normal modes) for no normal 
flow (fixed pressure) conditions. Of course, both scenarios reveal 
the presence of a continuous spectrum irrespective of the bound- 
ary conditions employed. It seems as though the inviscid pCf and 
SRI problems have properties and stability characteristics that are 
interchanged. 

4.4 Questions and a conjecture 

From a more physically motivated standpoint, we have experi- 
mented with these set of boundary conditions because they allow 
one to exert some comparative control between conditions. It is 
shown in Appendix lAlthat disturbances in the BE, subject to these 
boundary conditions, have a total disturbance energy, E, which 
evolves according to the exchange of energy that takes place be- 
tween the (Keplerian) shear and disturbance modes via the Orr- 
Mechanism and measured by the Reynolds Stress term, i.e. the 
RHS of <A3> . Conditions other than these would cause there to be 
some net work (positive or negative) to be performed on the layer 
during the ensuing course of the disturbances (Schmid & Henning- 
son, 2001). 

It is a puzzle, then, that in this QHSG limit there is an insta- 
bility in the case of no-normal flow conditions and none otherwise. 
Although this is merely a conjecture, is it possible that the SRI 
occurs because of the double reflecting boundary conditions? The 
instability shares many of the same properties of the acoustic insta- 
bility uncovered by Papaloizou & Pringle (1984,1985), otherwise 
known as the PP instability (Li et al. 2000). It is an instability of an 
acoustic disturbance in a domain like this with reflecting inner and 
outer walls in which there exists a critical layer, sometimes referred 
to as a corotation radius (Li et al., 2000). The waves grow in a reso- 
nant fashion because the reflecting walls, either one or both, allows 
for repeated passages of the wave across the critical layer which 
allows it to draw energy from the shear (Drury, 1985). It was found 
that the PP instability vanishes when the amplifying agent, usually 
the second reflecting boundary, is removed (Narayan, Goldreich & 
Goodman, 1987). 

Like the PP instability, the SRI as determined in this work are 
waves existing in a domain containing a critical point along with 
reflecting boundaries. When one of the boundaries no longer re- 
flects, there is no instability. Although these are neither compress- 
ible modes nor two-dimensional is it possible that the SRI arises 
in an analogous way due to the pathology that afflicts the PP insta- 
bility? This is an open question which should be clarified in future 
work. 

A clue towards this end might be found in the observation that 
there exists a second energy integral, as developed in Appendix D, 
involving a total energy expression T which says something inter- 
esting. The domain integral of the quantity T is conserved under 
no-normal radial flow conditions whereas it is not for the others. 
Perhaps the instability is related to this constraint placed on the dy- 
namics of the system? 

4.5 Flow two-dimensionalization and another conjecture 

We demonstrated in Section 2. 1 that there exists a conserved quan- 
tity (H) of the general linearized system of the BE that is adverted 
by the local basic shear. This quantity, which looks like a potential 
vorticity, is conserved independent of the QHSG asymptotic limit 
explored. Its analog is implied to exist in the for the LSB model 
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equations as discussed by Tevzadze et al. (2003). According to its 
definition, 1241 . 3 is composed of the vertical vorticity and a quan- 
tity representing buoyancy motions driven by density fluctuations. 
We also noted that in the limit where the buoyancy oscillations be- 
come very strong the term associated with it in 3 may become less 
important. In this circumstance it implies that the flow will take on 
a nearly two-dimensional character. 

In particular if the quantity iV 2 becomes large then, according 
to {19} one possible scaling between quantities in an initial value 
problem is to have O remain an order 1 quantity while the vertical 
velocity, w, be O (iV~ 2 ) . If all other quantities remain correspond- 
ingly order 1, that is to say if u, v, P, d x ,d y ,d z , dt ~ O (1), then 
to lowest order it would imply that the disturbances are dynamically 
in hydrostatic equilibrium and it would imply that the flow is nearly 
two dimensional conserving its vertical component of vorticity (cf. 

Barranco and Marcus (2005) demonstrated, in their shearing 
sheet simulations of a stratified fluid, the appearance of coherent 
vortical structures with vorticity vector pointing in the vertical di- 
rection. When they manifested themselves, the anticyclonic vor- 
tices appear near the vertical boundaries of the system, in other 
words, in that part of the atmosphere where the vertical component 
of gravity is greatest in magnitude. They also demonstrated the ro- 
bustness and persistence of these anticyclones by artificially remov- 
ing the vortex structure (after having developed) and replacing the 
flow field with noise. They show that the noisy spectrum quickly 
redeveloped into coherent anticyclone(s) much as it is known to do 
so in two dimensional shear flows (e.g. Umurhan & Regev, 2004). 
This fact is consistent with the implications suggested by the ad- 
verted conservation of the linear quantity 3. Of course, only a non- 
linear reformulation and reexamination of 3 can offer a more solid 
basis to any connection that there may exist here. 

Is it possible that it is a generic feature of stratified flow with a 
Couette shear profile and a vertically dependent Brunt- Vaisaila fre- 
quency (e.g. appropriate for a local representation of a circumstel- 
lar disk as here) to behave two-dimensionally in substantial parts 
of the atmosphere significantly away from the disk midplane, i.e. 
those regions dominated by a large Brunt- Vaisaila frequency? 
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APPENDIX A: AN INTEGRAL STATEMENT FOR THE 
BOUSSINESQ EQUATIONS 

It is instructive to consider integrals of the system since they can 
help guide one into deciding which boundary conditions to be used. 
We begin with by noticing that for the situations considered in here, 
the functional forms relating g and d z Tb are always constant multi- 
plicative factors of each other (see Sections 2.1 and 2.2). Therefore 
we take the ratio of these two quantities to be always a constant, 
that is, 

dzTt/g = constant, 

over the full spatial domain under consideration. With this assump- 
tion in hand one may (i) take the scalar product of 181101 and pt,u' 
, (ii) multiply lilt by 8ga p /d z Tt, and (iii) adding the results of (i) 
and (ii) together and making use of the incompressibility condition 
{7) to reveal 

(d t - qttoxdy) £ + v ■ V {£ + p) = 0, (Al) 

where 

= Pftu' 2 , gay Q 2 
2 d z T b 2 ■ 

Using condition Q once more, we may integrate <AH over the full 
spatial domain to find, 
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dE 
~~dt ' 

with 

E = 



Js 



ndS,- 



qu'v'dV 



(A2) 



d z T b 2 



dV, 



in which V and S is the volume and surface-boundary of the do- 
main in which n is the unit normal of the surface. The above result 
is true in general for both linear and nonlinear perturbations. We 
interpret the quantities in £ in the following way: the term p^u' 2 /2 
represents the kinetic energy in the disturbances while the term 
goL p 8 2 /2d z Tb represents the energy in thermal processes. By defi- 
nition £ is zero in steady state, while the steady pressure is constant, 
denoted by p. 

The global integral E, which we refer to as the total energy in 
disturbances, can change due to the influx of £ across the bound- 
aries, through work done upon the system externally as represented 
by the boundary flux term J s pu' ■ &dS, and finally due to the in- 
teraction with the background shear via the Reynolds stress term 
— J qu'v'dV (for a discussion of this see Schmid & Henningson, 
2002). 

The general evaluation of <A2> may proceed once boundary 
conditions are specified. As we have stated earlier, we will con- 
sider disturbances to be periodic in the y and z directions (sinuous 
or varicose for the latter). The radial boundary conditions and the 
motivation for their choices deserve some additional reflection. We 
remind the reader that one of the goals of this study is to assess 
whether or not the SRI is an intrinsic instability of the fluid and 
not some artifact of boundary conditions. One reasonable control 
is to require that there is neither work done on the system from 
outside nor there be any flux of energy across the bounding walls. 
This requirement requires that either the normal velocities are zero 
on either of the two walls or that (for linear disturbances only) the 
Lagrangian pressure perturbations are zero at the two bounding sur- 
faces. A mixture of these can also affect the same outcome. That is 
to say, for example, one could require that the normal velocity at 
one bounding surface is zero while the Lagrangian pressure per- 
turbation is zero at the other surface. Imposing these conditions 
therefore implies that E can change only due to the interactions of 
the perturbations directly with the shear, in other words, all such 
disturbances behave according to 



dE 
~dt 



(A3) 



In this sense, then, these solutions share some common property 
that allows for some comparison. 



APPENDIX B: DEVELOPMENT OF THE CONTINUOUS 
SPECTRUM MODE 

The discussion here largely follows the tack taken by Case (1960). 
Beginning with J35I we can say that 



(8t - Ftp 2 ) P = 0, 
is true for x 7^ x c where 



(Bl) 



qfloati 

Taking P^ 1 to denote the solution of IB II to the left and right 
(respectively) of x c , and, assuming zero pressure conditions at 
x = 0, 1 we have that 



A continuous spectrum mode: to = 0.0227 




Figure Bl. An example of the numerically generated continuous spectrum 
mode for fixed pressure conditions with a = 0.1, q = 3/2, On = 1. 
This particular mode has a frequency ui 0.023. Every seventh data point 
is plotted for the sake of visual clarity. The analytic representation of the 
eigenmode, i.e IB3I . is shown in the top graph for pressure only. In all cases 
the pressure eigenfunction is normalized to be 0.01 at x = 0.99. The jump 
in the horizontal (azimuthal) velocity is plainly evident. 



P = A sinh [k F x] . 



P + = A" 



sinh (1 - a;)] (B2) 



Then, to enforce continuity of the pressure at x = x c we see that 
sinh [k F (l — x c )] 



A 



A^ 



sinh [LiJ 



Once some normalization is specified, that is, a value of A is as- 
sumed, the solution is complete. For numerically generated solu- 
tions, we set P(x = 0.99) = 0.01. In summary, then, we have the 

(c) 

continuous mode pressure eigenfunction, P , is 

sinh[fc J? , (1 — x c )] 



P (c) = A + ■ 



sinh [kp x c 

sinh [k p (1 



■ sinh \k„x] < x < x c 



x)] 



(B3) 



X c < X < 1. 



in which A + = 0.01/ sinh [0.99 x k F ]. Note that it means that 
any value of U\ which satisfies the requirement < x c < 1 is 
an allowed solution. This is the nature of the continuous spectrum 
(Schmid & Henningson, 2001). Also, the mode associated with the 
continuous spectrum is not technically defined at x = x c . 



APPENDIX C: SUBTLETY IN THE K v OSCILLATING 
MODE SOLUTION 

A closer analysis of \55\ reveals that it is everywhere analytic and, 
therefore, an entire function along the real z axis, despite the pres- 
ence of a branch point at z — for the K v Bessel function. It 
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means, therefore, that in order to express the nature of this function 
as you one passes through z — one must be very careful about 
the relative phases that are incurred by crossing the z = point. To 
be more specific, let us define 



c = 



1/2 



and rewrite the solution 1551 according to its composition of Mod- 
ified Bessel Functions of the first kind, 



U(z) ~ C 5 



V) 



which, according to the series representation of l v would be 

= C f [c*T(3/4,C)-C*T(-3/4,C) 



where 



^ k\T(v- 

k=0 



1)' 



(Abramowitz & Stegun, 1972). The function T is symmetric with 
respect to the reflections f — £. However, the pre-f actors appear- 
ing in the above expressions imply that one must be very careful in 
interpreting the function as £ crosses over zero. This is best illus- 
trated by considering the behaviour of J55I for £ < 0. Expressed 
initially without restriction of £ J55I is, 

= C 3 T(3/4,C)-T(-3/4,C) 

Then if we restrict attention to £ < by defining the variable s = 
— £ and restricting attention to s > we see that the above becomes 

= - S 3 T(3/4, S )-T(-3/4, S ) 

-s*T(3/4, s) - s"'T(-3/4, s) 

-I (Is 2 ) -I (Is 2 ) 

-2T 3 (Is 2 ) +/C 3 



> 2 ) 



rewriting the above we see that the solution n(z) for z < be- 
comes as it is expressed in the text J56t . 



APPENDIX D: QHSG LINEARIZATION OF THE LSB 

We linearize <57I61> . It now means that p' and p' are the linearized 
density and pressure fluctuations. The resulting equations become 

(dt — qVL xd y )p + d x m u + d y m v + d z m w = 0, (Dl) 

= 2Qom v — d x p' , (D2) 

(d t — qCloxd v )m v + (2 - q)Q.om u = —d y p, (D3) 

= -d z p - p'g(z), (D4) 

(d t — qQoxd y )p b T,' + m w d z S b = 0, (D5) 

where for the sake of compact notation we have introduced the per- 
turbed momentum fluxes 



m u = p b u , 



m v = p b v , 



p b w 



and the perturbed entropy 



Operating on <D3> with d x and then making use of JO 1 1 reveals 



(dt — qSlaxdy) \d x m v — (2 — q)Sl p'~\ = (2 — q)Q,od z m w . 

This is followed by multiplying <D5> by fio (2 — g) /d z S b and then 
operating on the result with d z . This gives 



(d t - q n oX d y )^-n (2- q )^- 

Adding these two equations gives 

(dt — qClaxdy) 



(2 - q)Q d z m u 



d x m v - (2 - q)Q p + S~^o(2 - q)jr^~ 



= o. 



Given the definition of £', the hydrostatic relationship <D4t . and 
the radial geostrophic balance <D2> recovers J63L 



APPENDIX E: A SECOND INTEGRAL STATEMENT OF 
THE BOUSSINESQ EQUATIONS 

Following the steps in SectionlXlone may generate a second energy 
integral. The dynamical equations J7I1 1> describe the evolution of 
the both the disturbance velocities, i.e. the velocity fluctuations 
over and above the steady state Keplerian velocity, and the tem- 
perature fluctuations. We denote the total velocity of disturbances 
in the frame of the shearing box as U and given to be 



U = — qQoxy + u = {u, v — qflox, w'} 



(El) 



As such the governing equations of motion 181 1 1 i are more con- 
cisely written in vector form as 

dtU + U-VU = -4-Vp 

Pb 

+ 2f2 z x (U + qSloxy) + ±ga v 9z (E2) 



d t 9 + U ■ V9 = -wd z T b 
V-U = 0. 



(E3) 
(E4) 



As we have posited d z T b and g multiplicative factors of each other 
over the full spatial domain under consideration. With this assump- 
tion in hand one may (i) multiply JE2> by p b v , (ii) multiply IE3I 
by 9ga p /d z T b and (iii) adding the results of (i) and (ii) together to 
yield 



duF + U- V(.F + p) = 0, 



(E5) 



where 



d z T b 2 



n 2 2 
qilgX . 



With use of the incompressibility condition <E4t we may integrate 
<E5i over the full spatial domain to find, 



- = -J s (f + P )V.MS, 
with 



(E6) 



$ : 



TdV 



p b U 2 ga v 8 2 2 2 > ,, • 

— + dMY- qn ° x ^ dV > 



in which V and S are as they were defined before. We interpret 
the quantities in T in the following way: the term p b \J 2 /2 repre- 
sents the kinetic energy, the term — qV^x 2 is like a potential en- 
ergy and ga p 8 2 /2d z T b represents the energy in thermal processes. 
The global integral $ can change due to the influx of T across 
the dynamically undulating boundaries as well as through the work 
done upon the system externally as represented by the boundary 
flux term J s pU • ndS. 
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The point of this exercise is to note that only for no-normal 
flow boundary conditions does $ remain fixed for disturbances. 
The other conditions, like fixing the Lagrangian pressure, can cause 
$ to vary over the course of its evolution. This is because although 
p may be constant in steady state, the total quantity T is not con- 
stant in steady state where a simple inspection of its definition 
clearly reveals. By fixing only the Lagrangian pressure fluctuation, 
the otherwise moving boundary can allow T to seep in and out of 
the domain. 

Perhaps, then, the reason for the existence of the instability 
under no-normal flow boundary conditions arises because of this 
preserved property of the disturbances. The reflection property of 
the boundaries perhaps traps energy in a way that causes growth 
to be encouraged. In this case, the energy of the disturbances must 
come from the energy contained in the background shear state and 
because there is an overall trapping of the energy, a runaway ex- 
traction processes takes place - ironically, leaving the total energy 
budget , $, fixed over the course of the evolution. 



